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Abstract 

The flow of a simple glass forming system (a 80:20 binary Lennard-Jones mixture) through a 
planar channel is studied via molecular dynamics simulations. The flow is driven by an external 
body force similar to gravity. Previous studies show that the model exhibits both a static [Varnik 
et al. J. Chem. Phys. 120, 2788 (2004)] and a dynamic [F. Varnik and O. Henrich Phys. Rev. 
B 73, 174209 (2006)] yield stress in the glassy phase. These observations are corroborated by 
the present work, where we investigate how the presence of a yield stress may affect the system 
behavior in a Poiseuille-type flow geometry. In particular, we observe a blunted velocity profile 
across the channel: A relatively wide region in the channel center flows with a constant velocity 
(zero shear rate) followed by a non linear change of the shear rate as the walls are approached. The 
observed velocity gradients are compared to those obtained from the knowledge of the shear stress 
across the channel and the flow-curves (stress versus shear rate) , the latter being determined in our 
previous simulations of homogeneous shear flow. Furthermore, using the value of the (dynamic) 
yield stress known from previous simulations, we estimate the threshold body force for a complete 
arrest of the flow. Indeed, a blockage is observed as the imposed force falls below this threshold 
value. Small but finite shear rates are observed at stresses above the dynamic but below the static 
yield stress. We discuss the possible role of the stick-slip like motion for this observation. 

PACS numbers: 64.70.Pf,05.70.bn,83.60.Df,83.60.Fg 



I. INTRODUCTION 

The so called soft glassy materials [l, 2, 3] exhibit a rich variety of interesting rheo logical 
phenomena. When compared to a Newtonian liquid (defined as a liquid whose viscosity, 77, 
does not depend on shear rate, 7) the viscosity of a soft glassy system shows pronounced 
dependence on imposed shear rate. To be specific, let us consider disordered colloidal sus- 
pensions typical example. It is well known that the shear viscosity of dense 
colloidal dispersions decreases with increasing 7 (shear thinning) if one focuses the attention 
on low shear rates. In the limit of high shear rates, on the other hand, the shear viscosity 
starts to increase with 7 (shear thickenning). While shear thinning is commonly attributed 
to the competition between the time scale imposed by the external flow and the time scale 
of inherent structural relaxation, shear thickenning phenomenon is rather understood as 
originating from hydrodynamic effects |?| (whose contribution to the stress is negligible at 
low shear rates and high concentrations but increases considerably at high 7 jsj]). 

In this paper, we study a model system whose rheological properties can be rationalized 
without taking into account hydrodynamic effects [9j. Previous studies of the model showed 
that it exhibits both a static and a dynamic yield stress. The main difference between the 
static and the dynamic yield stress lies upon the imposed quantity. While the static yield 
stress, (j^ tatlc ; is measured in experiments upon imposed stress, the dynamic yield stress, 
<7y ynamic , is measured in experiments upon imposed shear. This is so because a soft glassy 
material does not develop the same rheological response in the both mentioned cases. 

To see this, let us consider a planar Couette cell. If a lateral force per unit area (=stress) 
is imposed to one of the walls of the Couette cell, the system may resist to the imposed stress 
if the latter is below some threshold value (which generally depends on temperature and 
density, say). For stresses higher than this threshold value, on the other hand, the system is 



10] . The static yield 



'liquidized' and flows with a linear velocity profile across the channel 
stress thus characterizes the response of an initially non-driven amorphous solid. 

If instead of the stress an average shear rate is imposed (by, e.g., moving one of the 
walls with a constant velocity), the occurrence of a flow is unavoidable by construction. 
However, the flow profile need not be linear in this case and may exhibit strong spatial 
heterogeneity. In particular, a two-phase scenario may occur: A region of zero shear ('solid- 
like' or 'jammed') coexisting with a sheared ('liquid-like') region 11J. Interestingly, as the 
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wall velocity approaches zero, the shear stress across the system does not converge to zero 



(as would be the case for a Newtonian fluid) but seems to saturate at a finite value [9l.llll. Il2|. 
This limit of the shear stress for vanishing shear rate is usually defined as dynamic yield 
stress. It follows from the above description that the dynamic yield stress characterizes the 
response of a shear molten amorphous solid. 

It is noteworthy that the static yield stress is found to be higher than its dynamic coun- 
terpart [lol ]. a situation reminiscent of the difference between static and dynamic friction 
(without any claim for a strict analogy). For the present model and at a temperature and 



density of T = 0.2 and p = 1.2 our previous studies yield ^ ynamic « 0.5 Q and <xf atic « 0.6 
jl(3 | (all quantities are given in Lennard- Jones (LJ) units; see below). 

In a glass forming system, the static yield stress does in general depend on the system 
history and the way the stress is imposed. In particular, it depends on the waiting time (the 
time elapsed between the temperature/density quench and the beginning of stress ramp) as 
well as on the rate with which the stress is increased. In the present simulations, however, 
no stress ramp is applied. Rather, after a waiting time of t w , we instantaneously switch on a 
constant external force field exerting a body force of F e on each particle (see also section ILT1) . 
Furthermore, we focus on sufficiently large waiting times so that, within the time window 
accessible to our simulations, flow profiles are hardly affected by a dependence of the static 
yield stress upon t w (see section [TTT| for further discussion of this issue). 

In contrast to the static yield stress, the above definition of a dynamic yield stress allows 
one to avoid these complications. This follows from the fact that < jy yn3,mic reflects the sys- 
tem response to a small but finite steady state shear, where time translation invariance is 
recovered. Taking the limit of vanishing shear rate does not affect this property. 

Obviously, the existence of a dynamic yield stress presupposes the existence of a plateau 
in the flow curve (shear stress versus shear rate) in the limit of low shear rates. However, 
although our previous studies support the existence of such a plateau in the case of the 
present binary LJ model, the issue of a dynamic yield stress in soft glassy materials still 

n 

remains controversial (see e.g. [13J for recent experimental studies of this topic in colloidal 
dispersions). Therefore, it is worth to check whether our model system also exhibits other 
features which follow from the existence of a yield stress. As will be shown in this report, 
the answer to this question is affirmative. In particular, we observe non trivial behavior 
such as profile blunting and flow blockage in a Poiseuille-type flow geometry, features which 
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can consistently be described assuming the existence of a dynamic yield stress. 

The paper is organized as follows. After an introduction of the model and the simulation 
method in the next section, the effect of yield stress on the behavior of the system in a planar 
channel flow driven by an external body force will be investigated. In particular, it will be 
shown that the velocity profile exhibits salient features of a two-phase system: A 'solid-like' 
central part with zero velocity gradient and two lateral 'liquid-like' sections between the 
channel center and the walls. A consequence of this property on the dependence of the 
mass flow rate upon the imposed force is worked out and compared to the case of the same 
system in the normal liquid state, where it behaves like a Newtonian liquid with a shear 
independent viscosity. A summary compiles our results. 



II. A BINARY LENNARD-JONES MIXTURE 



In order to address the above mentioned issue, we study via molecular dynamics simula- 
tions a generic glass forming system, consisting of a 80:20 binary mixture of Lennard- Jones 
particles (whose types we call A and B) at a total density of p = Pa + Pb — 1-2. 

A and B particles interact via U^{r) = 4e Q/3 [((i C j / g/r) 12 — (d a p/r) e ], with = A,B, 
6ab = 1-5caA) eBB = 0.5eAA, g?ab = 0.8<iAA, g?bb = 0.88<iAA and mB = mA- In order to enhance 
computational efficiency, the potential was truncated at twice the minimum position of the 
LJ potential, r C)a p = 2.2A5d a p. 

The truncation of the LJ potential introduces a discontinuity in the force field, which 



could be corrected via a smoothing procedure 



14|. However, the present model with the 



truncated version of the LJ potential has extensively been studied in the literature and has 
become a benchmark model for the studies of glassy systems. Therefore, and for the purpose 
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of comparison with previous studies [lOj, [Hi, Ilia, [la, [17|, lis, llSj, we keep the model as it is. 
Note also that the use of a truncated LJ potential is not a priori a disadvantage, since 
we do not seek a comparison with analytic studies of this specific system. Rather, we are 
interested in generic features of a glass forming model system for which the present binary 
LJ mixture has indeed become a prototypical example. 

The parameters 6aa, g?aa and ttt-a define the units of energy, length and mass. All other 
quantities reported in this paper are expressed as a combination of these units. The unit of 
time, for example, is given by tlj = ^aav^a7^aa and that of stress by caa/^aa- Equations 
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of motion are integrated using a discrete time step of dt = 0.005. In order to test numerical 
accuracy, we also performed simulation runs using a smaller time step of dt = 0.001. Since 
no deviations were found between the results obtained for dt — 0.005 and dt — 0.001, we chose 
the larger time step for all the subsequent simulations. 

The present model has first been introduced by Kob and Andersen in the context of the 
dynamics of supercooled liquids 
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171 ] , who showed that it was suitable for an ana 



ysis 



of many aspects of the mode coupling theory of the glass transition (MCT) [2d, l2l|, |22|. 
In particular, at a total density of p = 1.2, equilibrium studies of the model showed that 
the growth of the structural relaxation times at low temperatures could be approximately 
described by a power law as predicted by the ideal MCT, r rc i ax oc (T — T C )~ 7MCT . Here, 
T c = 0.435 is the mode coupling critical temperature of the model and 7mct is the critical 
exponent. For the present binary Lennard- Jones system, numerical solution of ideal MCT 



equations yields a value o: 
mixture of soft spheres 



7mct 



i.5 Q. 



A similar value is also obtained for a binary 



The model has also been studied in the context of the rheology of disordered systems 
Q, H, 11, [l9, l24j . In these studies, aspects such as flow heterogeneity [11] . structural 




FIG. 1: (Color online) A snapshot of the simulation box. The system consists of 80 percent of 
particles of type A (blue) and 20% of type B (red). Two atomistic walls (yellowish colors) confine 
the system along the z-direction (to the left and right of the image). Periodic boundary conditions 
are applied in the remaining x and y directions. The system size is L x x L y x L z = 30 x 30 x 86. 
It contains 92880 particles in total. The center of the coordinate system (r = 0) is placed in the 
middle of the simulation box. Thus, x G [—L x /2 L x /2], y £ [—L y /2 L y /2] and z € [— L z /2 L z /2]. 
x is the flow direction, y the neutral (vorticity) direction and z the direction of the velocity gradient. 
All lengths are in LJ units. 
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191 ] and the existence of a static [lOj and a dynamic [9J yield 
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relaxation under external drive 
stress were addressed. 

An interesting consequence of the presence of a yield stress on the flow behavior of the 
system is presented here. For this purpose, we simulate a Poiseuille-type flow in a planar 
channel. The study of such a situation is interesting since the stress in a Poiseuille-type flow 
is zero in the channel center and increases linearly with the distance from it. As will be shown 



below, this is a consequence of momentum balance equation [25|, [26j and thus independent of 
the specific flow profile formed across the channel. If the fluid under consideration exhibits 
a finite yield stress, one may expect that the fluid portion in a certain region around the 
channel center (where the stress is below yield stress) should behave like a solid body while 
it should flow like a liquid further away from this region. 

We first equilibrate a system of size L x x L y x L z = 30 x 30 x 86 (containing 92880 
particles) at a temperature of T = 0.45(> T c ). The temperature is then set to T = 0.2 
(< T c ) corresponding to a glassy phase. Note that, the velocity distribution adapts itself 
very fast (within a time of one LJ unit) to the Maxwell distribution corresponding to the 
new temperature. Particle configurations, on the other hand, keep the memory of the old 
(high) temperature for far larger times (see e.g. [10] for a short discussion of aging in the 
case of the present model system). 

However, as time proceeds, particles gradually rearrange and the system moves towards 
that part of the configurational space which corresponds to the new (low) temperature. This 
process is fast in the beginning, where thermodynamic driving forces are the largest, but 
slows down with time. In fact, as a characteristic feature of a glassy phase, the system never 
reaches thermal equilibrium. Rather, it keeps evolving towards it endlessly (aging). 

As a result of this aging process, time translation invariance is violated and those system 
properties which would be independent on the measurement time in an equilibrium state 
(such as structural relaxation time, diffusion coefficient, static structure function, etc.) show 
a dependence on the time elapsed between the temperature quench and the beginning of the 
measurement. Furthermore, quantities computed as time averages also show a dependence 
on the duration of the measurement, thus reflecting the fact that, in an aging system, 
ensemble and time averages are no longer equivalent. 

In particular, in a glassy state, inherent system dynamics slows down upon aging and 
structural relaxation times grow (ideally) endlessly, eventually exceeding any other time 
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scale in the problem (such as the time scale imposed by external shear). In this interesting 
limit, the system no longer behaves like a liquid but rather exhibits solid-like properties, as 
exemplified by the presence of a finite static yield stress. Since we wish to concentrate on 
this late time behavior, we first run simulations in the quiescent state for a sufficiently large 
waiting time of t w = 10 4 LJ units before imposing an external force. 

A time of t w = 10 4 LJ units in our simulations is sufficiently large in the following sense. 



First, it is large so that the system exhibits a finite, measurable static yield stress [10]. 
Second, it is large enough so that the increase of the static yield stress upon further aging 
is slow and does not lead to a qualitative change in the flow behavior (see section II I II for a 
more detailed discussion of aging effects on the flow behavior). 

After this initial period of time, two solid walls are introduced parallel to the xy-plane 
by immobilizing all particles whose z-coordinate satisfies \z\ > 40 (walls of three particle 
diameter thickness). A flow is then imposed along the x-axis by applying on each particle 
a constant force, F e . This gives rise to a force density of / = pF c . The force on a particle 
is the sum of the interaction forces arising from the Lennard- Jones potential and F e (recall 
that F e = for t < t w ) . 

For our geometry with planar solid walls, it follows that (u-'V)u = 0, where u = (u(z), 0, 0) 
is the streaming velocity (recall that x is the flow direction, y the neutral (vorticity) direction 
and z the direction of the velocity gradient). In the steady state and in the absence of a 
pressure gradient (Vp = 0), the momentum continuity equation thus reduces to da/dz = 
f = pF c , which yields a(z) = pF c z, where we used the symmetry of the shear stress with 
respect to the rn/-plane (a(z = 0) = 0). 

The external force does work on the system. This work is transformed into heat via 
viscous dissipation. In the absence of a thermostat, this would lead to a continuous increase 
of temperature with time. In order to keep the system temperature at a prescribed value, the 
viscous heat must be removed. For this purpose, we divide the system into parallel layers of 
thickness dz — 1 and rescale (once every 10 integration steps) the ^/-component of the particle 
velocities within the layer, so as to impose the desired temperature T. More precisely, we 
first compute the local kinetic energy per particle e^ Q = l/N(z)Y^ m i v yi within a layer 
centered at z. Here, is the mass of a particle, N(z) the number of particles in the layer 
and the sum runs over the particles in the layer only. A scale factor, s, is then determined 
via the requirement that the new velocities sv y ^ satisfy 1/N(z) Yl m i{ sv yi) 2 = k~sT. This 
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gives s = {k-oT /e^) 1 / 2 . Finally, the new velocities are computed via multiplication of v y ^ 
with s. 

Note that such a local treatment is necessary to keep a homogeneous temperature profile 
when the velocity profile is not linear. This is so because in this case the shear rate, 7, 
and hence the rate of heat production, 7/7 2 may vary within the channel giving rise to a 
temperature gradient if only the average temperature in the channel is adjusted. 

However, despite this local temperature control, the heat production close to the walls 
(where the shear rate is very high) is so fast that a temperature increase in this part of the 
channel can not be fully avoided (recall that viscous heat scales with 7 2 ). The magnitude of 
the excess temperature strongly depends on the applied force per particle and is practically 
negligible for F c < 0.02 (see also the discussion of Fig. [6]). 

III. RESULTS 

Let us first examine possible effects of aging and flow time on the velocity profile. For 
this purpose, we prepare the system as described above for various choices of the waiting 
time t w . Recall that t w is the time interval between temperature quench and the time at 
which the external force is switched on. At this instant, we set the clock to zero, t = thus 
marking the onset of body force. 

While the system would gradually 'solidify' in the quiescent state, the external force tries 
to induce a flow and thus tends to 'fluidize' (rejuvenate) the system. The extent and the 
rate of this rejuvenation does, however, strongly depends on the magnitude of the stress 
formed across the channel. As mentioned above, the stress in the present Poiseuille-type 
geometry is a linear function of the distance from the channel center: a(z) = pF e z. For a 
given external force, the shear stress thus increases when approaching the walls. 

As a consequence, the system flows first in the proximity of the walls, while the center 
of the channel behaves rather like a solid body. The effect of aging now reflects itself in the 
onset of the flow. The larger t w the slower the initial flow behavior. This feature is nicely 
born out in the left panel of Fig. [21 As a survey of the data corresponding to t = 10 2 reveals, 
the system hardly flows for t w = 10 3 and t w = 10 4 (deforming rather like an elastic solid) 
while for t w = 10 2 it 'liquidizes' in the proximity of the walls. 
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The data shown in the left panel of Fig. [2] also suggest that, for larger waiting times, the 
external force must be applied during a longer period of time (larger t in the figure) in order 
to remove memory effects. This is seen by a closer look at the data corresponding to t — 10 3 . 
Here, the curves belonging to t w = 10 2 and £ w = 10 3 are practically identical, while that of 
t w = 10 4 shows significantly slower deformation behavior. However, applying the external 
force during a time of t = 5 x 10 3 , aging effects disappear also in the case of t w = 10 4 . 

For a given distance from the channel center, the shear stress increases (and hence the 
time scale imposed by the external force decreases) with increasing the magnitude of the 
applied force. This is an important aspect since memory effects (such as effect of aging) decay 



within a typical structural relaxation time, the 
imposed time scale in glassy systems (see e.g. 



atter being of the order of the externally 
for a detailed discussion of this issue). In 



fact, aging effects are expected to be observable as long as this externally imposed time scale 
is significantly larger than the waiting time. Since this time scale reduces upon increasing 
F e , aging effects are expected to also disappear faster. This property is nicely born out in 
the right panel of Fig. El where data similar to the left panel are shown for the choice of a 
larger external force per particle, F e = 0.025. As seen from this panel, already at a time of 
t = 10 2 , the flow behavior of the system is practically independent of t w , even for a waiting 
time of t w = 10 4 > t = 10 2 . 

Unless otherwise stated, the waiting time is t w = 10 4 for all the results presented below. 
After the onset of external force, we wait another period of time (of duration t = 5 x 10 3 ) 
before starting the data collection. This helps to reduce the above discussed transient effects 
related to a competition of aging and shear. Note, however, that in a Poiseuille-type flow 
of a glassy system, transient effects will always be present to some extent. This is closely 
related to the fact that the shear rate approaches zero as one moves from the walls toward 
the center of the channel, thus giving rise to progressively large relaxation times. As a 
consequence, the time necessary to establish steady state ideally diverges in the center of 
the channel. This behavior is enhanced in a yield stress fluid, where the zero-shear zone has 
a finite extension comprising a region around the channel center where the local shear stress 
is below the fluid's yield stress. The 'jammed' region corresponds to this part of the system. 

In the case of a Newtonian liquid, the above mentioned approach to drive the flow would 
give rise to a parabolic velocity profile of the form u(z) = pF c (L 2 z / '4 — z 2 )/(2r]). However, as 
Figs. El and [3] clearly demonstrate, the velocity profile in the case of the present model in the 
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transverse coordinate, z transverse coordinate, z 



FIG. 2: (Color online) Competitive effects of aging and shear rejuvenation on the velocity profile. 
Left: Velocity profile for waiting times of t w = 10 2 , 10 3 and 10 4 (LJ units) and various measurement 
times, t, as indicated. Here, i w is the time elapsed between the temperature quench and the onset 
of forcing (F e = 0.0175). Data collection starts at times t = 10 2 , 10 3 and 5 x 10 3 after switching 
on the external field. In order to better reveal the effects of aging, the duration of a measurement 
is limited to dt = 10 2 , 10 3 and 5 x 10 3 , respectively. Right: Similar plot as in the left panel, but 
for a higher applied force of F c = 0.025. Here, data collection starts at times t = 10 2 and 10 3 after 
switching on the external field and spans over dt = 10 2 and 10 3 , respectively. Due to stronger 
forcing, effect of shear rejuvenation shows up at shorter times compared to the left panel. Note 
that the system response is independent of t w even at t = 10 2 (the observed difference between the 
curves is statistical, otherwise a smaller deformation would be observed at t w = 10 4 compared to 
t w = 10 3 ). All quantities are in LJ units. 

glassy phase exhibits a quite different behavior: In the central region, the velocity profile is 
flat with a zero gradient while it gradually departs from this constant behavior (shear rate 
becoming non zero) beyond this central part. 

It is interesting to note that similar (blunted) shapes of the velocity profile are also 
observed in pressure driven flows of both neutrally buoyant suspensions of spheres (with 



a size of the order of 1mm) |27l l28lj as well as red blood cells (biconcave disks of 2fim 



thickness and 8/im diameter) |29l. l30j| . In these cases, however, the profile blunting is usually 
accompanied by a migration of particles from the wall region towards the center of the 
channel ('wall migration') a phenomenon, which is absent in the case of present studies (see 
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the right panel of Fig. [3]). 

Other examples of blunted velocity profiles occur in systems with a nematic order pa- 
rameter (see e.g. 3J| and references therein). In these systems, the non-trivial rheological 
response is closely related to a variation of the system structure accompanied by a change in 
the nematic order parameter. However, such an order parameter-related structural change 
is absent in the case of our glass forming model. 

Nevertheless, a qualitative similarity to the case of present simulations may be found 
when red blood cells are concerned. This similarity rests upon the fact that profile bunting 
in red blood cells occurs only if a certain amount of aggregation among red blood cells is 
present (see e.g. Fig. 7 in reference 30J). The aggregation enhances the solid character of 
the suspension and leads to a higher yield stress, similar to a reduction of temperature in 
our model. 

The width of the "jammed" region can be estimated from a knowledge of the yield stress 
in the system via o{z = W/2) = pF c W/2 = a y which gives W = 2cr y /(pF c ). The question 
arises whether the static or the dynamic yield stress should be used for an estimate of W. 
Results of our simulations suggest that the dynamic yield stress yields a better estimate of 
the width of the solid-like region in the channel center. This is exemplified in Fig. [31 In 
this figure, two vertical dashed lines mark the bounds for the solid-like region estimated via 
dynamic yield stress whereas the bounds denoted by short vertical solid lines are obtained 
using the static yield stress. As a survey of the velocity gradient reveals, the use of cTy tatlc 
overestimates the width of the solid-like region. 

It is interesting to check the origin of the observed finite shear rates at stresses above the 
dynamic but below the static yield stress. For this purpose, we performed a series of long 
simulation runs of a smaller system (L x X Ly X L z = 10 x 10 x 40) at constant imposed 
stress for a temperature of T — 0.2 (far below T c ). All simulations started after an initial 
aging of the system during a time of t w = 10 4 . Interestingly, our data reveal the presence of 
a stick-slip like plastic deformation. This occurs not only at stresses between <jdynamic anc j 
^.static a i so ^ s t resses below <jdy namic . 

This behavior is illustrated in the left panel of Fig. HI The panel shows the center 
of mass position of the whole fluid versus time for some selected values of the imposed 
stress ranging from o = 0.46 (below the dynamic yield stress) up to the static yield stress 
(er = (j^ tatlc = 0.6). The corresponding center of mass velocity is depicted in the right panel 
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FIG. 3: (Color online) Left: Flow through a three dimensional planar channel of the model fluid 
studied in this paper. The flow is generated by imposing a constant body force of F e = 0.025 on 
each particle. At a temperature of T = 0.2 and a total density of p = 1.2, the model exhibits 
a stress plateau at low shear rates which, for the present purpose, plays the role of a dynamic 



yield stress, a 



dynamic 

y 



0.5 



One can easily show that the stress in the channel behaves as 



32l | . Obviously, 



a = pF c z, where z is the transverse coordinate (z = being the channel center) 
for \z\ < 16.67(~ cr y /(pF c )), the stress in the channel is below the dynamic yield stress. One thus 
expects the system to behave as a solid body for \z\ < 16.67: Either it should be at rest or move 
with a constant velocity (zero velocity gradient). Indeed, an inspection of the velocity profile, u(z), 
and its derivative (rescaled to fit into the figure) confirms this expectation (the region delimited 
by the two vertical dashed lines). For \z\ > 16.67, on the other hand, the local shear stress exceeds 



a 



dynamic 

y 



0.5 leading to a liquid like behavior. The system flows with a shear rate which non- 



linearly increases upon increasing stress. The short vertical lines show the limits of the expected 



'jammed' region using the static yield stress, <jS tatlc ~ q.6 (see also the text) 



2Jj. Right: Profiles 



of the system density and the component of the pressure perpendicular to the walls. Both the 
density and the normal pressure are practically constant across the channel. In particular, there 
is no wall migration in the case of present simulations. This strongly suggests that the blunting of 
the velocity profile observed in the left panel is related to the presence of a yield stress which leads 
to the mentioned two-phase behavior. All quantities are in reduced (LJ) units. 
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of the same figure. As seen from this panel, for stresses below a = 0.52, the contribution 
of the stick-slip like motion to the flow velocity drops by roughly an order of magnitude. 
It is noteworthy that a qualitatively similar stick-slip like dynamics has also been observed 
in our previous simulations under imposed shear llj as contrasted to the present case of 
imposed stress. 

Using the data shown in Fig. HI we can estimate the contribution of this intermittent 
motion to the overall shear rate via the relation 7 = V cm /L z . Using L z = 40, this gives 
*1 = 5.9 x 10- 4 , 7.3 x 10~ 5 , 6.1 x 10" 5 , 1.1 x 10~ 5 , 1.1 x 10~ 5 and 2.8 x 10~ 7 for a = 
0.60, 0.56, 0.52, 0.50, 0.48 and 0.46 respectively. The finite shear rate observed in the case 
of stresses between a = 0.5 and a = 0.6 is therefore closely related to the onset of significant 
stick-slip motion at these stresses. 

One interesting consequence of the presence of a yield stress is that, for a given driving 
force (pressure gradient) a flow blockage may occur if the channel width is too small to 
allow the formation of stresses above the system's yield stress. This 'critical' channel width 
is simply estimated via L z = 2a y /(pF c ). Obviously, a similar situation would also occur at 
constant channel width via decreasing the applied body force below F c = 2a y /(pL z ). Let 
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FIG. 4: (Color online) Left: X-position of the center of mass of the fluid versus time for some 
selected stresses. The results shown here correspond to simulations under constant imposed stress 
of the present binary LJ model at a temperature of T = 0.2 (below T c ). The dimensions of the 
simulation box are L x x L y x L z = 10 x 10 x 40. Right: The same data as in the left panel, now 
divided by time. Note the presence of a flow with a time-independent velocity at a = 0.6. All 
quantities are in LJ units. 
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us estimate this 'threshold' F e . At the temperature and density studied here (T = 0.2 and 
p = 1.2) the system exhibits a dynamic yield stress of Uy yna,mic ~ 0.5. Using this value as 
the yield stress along with L z = 80 we obtain F e pa 0.01 for the minimum force per particle 
required to induce a flow in the system. 

This aspect is illustrated in the left panel of Fig. [5J where velocity profiles are shown 
for various choices of F e . As expected, the width of the 'jammed' region increases with 
decreasing F e . For F c = 0.017, for example, there remains a narrow liquid-like region close 
to the walls while the rest of the system is in a 'jammed' state. This observation is made 
more quantitative in the right panel of the same figure. For this purpose we determine for 
each F e the width of a region with a shear rate larger than a small but finite value. We find 
that the choice |7| = 10 -3 is a good choice for an accurate determination of the position of 
the 'interface' between the liquid-like and the solid-like regions. As can be seen from the 
right panel of Fig. [5], results of this analysis obey well the expected relation W = 2a y / (pF e ) if 
cr y = ^dynamic ~ o.5 is used. The use of static yield stress (<7y tatlc ~ 0.6) would overemphasize 
the size of the solid-like part of the channel. 

The left panel of Fig. [6] illustrates the velocity gradients, du(z)/dz, for exactly the same 
values of the external force per particle as shown in the left panel of Fig. [5j Here, velocity 
gradients are compared with profiles of shear rates, 7(17(2:)), estimated from a knowledge of 
the local stress and the flow curve upon homogeneous shear: For each z, the corresponding 
stress is first evaluated via a = pF c z. In order to estimate the shear rate corresponding to 
this shear stress, we use the homogeneous flow curves (7 — a) determined in our previous 
simulations [9] . In these simulations, an algorithm capable of ensuring a constant shear rate 
across the system was used (see e.g. [24j for details of the simulation method). 

As the inset of the left panel of Fig. [6] demonstrates, du(z)/dz and 7(17(2;)) agree well 
within a certain region in the channel, whose extension increases upon decreasing the exter- 
nally imposed force, F e . This trend is probably related to the variation of the temperature 
across the system. 

Indeed, the right panel of Fig. [6] shows that the system temperature is not constant overall 
in the channel but slightly deviates from the prescribed value in a region between the channel 
center and the walls. The extent of the central region with a constant temperature increases 
at lower F e while at the same time the magnitude of the excess temperature with respect to 
the prescribed value decreases. Noting that significant temperature excesses occur only as 
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the shear stress reaches values comparable to twice the dynamic yield stress of the model, 
it is not surprising that, at such high stresses, even a local thermostat is not able to ensure 
a constant temperature profile across the system. 

In the above, the knowledge of the homogeneous flow curve is used in order to obtain 
an independent estimate of the local shear rate across the channel. Similarly, one can use 
the fact that the shear stress is well known in the channel, a(z) = pF e z, along with the 
knowledge of the velocity gradient in order to obtain the flow curves, i.e. shear stress as a 
function of shear rate (velocity gradient). For this purpose, we plot in Fig.[7|for each value of 




transverse coordinate, z inverse force per particle, 1/F 



FIG. 5: (Color online) Left: Velocity profiles (LJ units) obtained for various choices of the external 
force per particle as indicated. At the temperature and density studied here (T = 0.2 and p = 1.2) 
the system exhibits a dynamic yield stress of a y ~ 0.5. Using this value as the threshold stress for 
the onset of shear flow along with the relation cr(z) = pF e z for the local stress in the channel, the 
expected width of a central region with a solid-like behavior (flat velocity profile) can be estimated 
via W = 2a y /(pF e ). The short vertical dashed lines help to illustrate this central region. Beyond 
this part of the channel, a liquid-like behavior is expected. As the external body force reaches a 
value of F c ~ 0.01, the 'jammed' region comprises the whole channel. Right: The width of the 
solid-like region is estimated approximately as the width of a region with a shear rate larger than 
7 = 10~ 3 (see the left panel of Fig. [H] for an inspection of shear rates). The solid line gives the 
above mentioned expected linear dependence on 1/F C . Deviations at high F c are probably related 
to undesired viscous heating which leads to an enhanced softening of the solid region (see also a 
discussion of the velocity gradients). All quantities are given in LJ units. 
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z the corresponding values of a(z) versus 7(2). Due to the symmetry around the mid-plane 
of the channel, the above procedure would yield identical flow curves using the data from 
the left and right halves of the channel provided that no statistical uncertainty is present. 
In reality, however, there is a finite statistical scatter. In order to illustrate this fact, we do 
not average the results but depict the individual flow curves obtained from the analysis of 
each half of the channel. The so obtained a — 7 curves are then compared to the result of 
simulations at homogeneous shear |9j . The comparison is restricted to shear stresses, where 
the relative deviations between the local temperature and the prescribed one are below one 
percent. With this restriction, the flow curve obtained from the present simulations agree 
well with that of homogeneous shear. 

Finally, Fig. [S] illustrates how the mass flow rate, Q = 2pL y j^ z ^ 2 u(z)dz, varies with 
applied force per particle. As a survey of the velocity profiles (Fig. [5]) already suggests, the 
mass flow rate rapidly decreases as the 'critical' force per particle F e m 0.01 is approached 
(note the logarithmic scale for the y-axis in Fig. [8]). In order to highlight the effect of a 
finite yield stress, the data corresponding to the same model in the normal liquid state 
(where the yield stress is identically zero) are also shown. For the normal liquid state, it is 
straightforward to show that Q oc F e . This relation is nicely born out by by the simulated 
data. We are, however, not aware of an analytic expression for the mass flow rate in a yield 
stress fluid. Therefore, we fit the data to the simplest polynomial in F c — 2cr y / (pL z ) which 
describes the simulated data best. 



IV. SUMMARY 



In this paper, we report on the flow behavior of a simple glass forming system. The model 
consists of a 80:20 binary Lennard- Jones mixture first introduced by Kob and Andersen in 



the context of the dynamics of supercooled systems 15 
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17J] - It is we 



capability to form a disordered solid at low temperatures/high densities 33]. 



1 known for its 



Previous studies of the rheological response of the model suggest the existence of a finite 
yield stress, <r y , in the glassy phase. In particular, at a temperature of T = 0.2 (deep in 
the glassy phase), a finite static yield stress of <j^ tatlc « 0.6 was found via simulations upon 
imposed stress lOj. More recent studies of the same model under imposed shear, on the 
other hand, showed the existence of a stress plateau in the low shear rate limit of the flow 
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FIG. 6: (Color online) Left: Shear rate, j(z) = du(z)/dz, for exactly the same choices of F c 
as in the left panel of Fig. [5j As expected, the width of the central solid-like region (7 ~ 0) 
increases upon decreasing the external force per particle and eventually reaches the width of the 
channel at F e ~ 0.01 (see also the left panel of Fig. [5]). The dotted dashed lines to the right of 
the panel represent velocity gradients computed using the flow curve (a — 7) of the same model 
under homogeneous shear [9|] and the fact that u{z) = pF e z (i.e., by plotting for each z the shear 
rate A f{a{z))). As seen from the inset, for z < 20 the thus obtained 7(2) agrees well with the 
observed velocity gradient. However, significant deviations occur close to the walls (main figure). 
These deviations as well as the maxima in \du(z)/dz\ are probably related to a local increase of 
temperature. A vertical dashed line at z = 38 marks the approximate position of these maxima. 
As a survey of the right panel reveals, temperature profiles also exhibit maxima at approximately 
the same transverse position. Right: Temperature profile across the channel for various choices 
of the external force per particle as indicated. The temperature profiles develop maxima in the 
proximity of the walls. The magnitude of the excess temperature with respect to the prescribed 
value (T = 0.2 in the present case) increases at higher F e . In addition to a stress-induced decrease 
of shear viscosity (shear thinning) this temperature increase enhances the decrease of the local 
viscosity further. As a result, the local shear rate, 7(2), increases faster than would be expected on 
the basis of shear thinning alone. As in the case of the left panel, a vertical dashed line at z = 38 
marks the approximate position of temperature maxima. All quantities are given in LJ units. 

curve (stress versus shear rate) in the glassy phase Q . As the present work also underlines, 
this stress plateau paly the role of a (dynamic) yield stress. 
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FIG. 7: (Color online) Shear stress versus shear rate. Filled symbols correspond to the data 
obtained from an analysis of the relation between the local shear rate and the local shear stress 
across the channel in the present simulations (F e = 0.0167). Squares correspond to the data 
extracted from the left part of the channel and triangles to the right half of the simulation box. 
Connected circles reproduce the flow curve obtained from previous simulations of a homogeneous 
shear for the same model and at the same density and temperature as studied here (T = 0.2, 
p = 1.2) A horizontal dashed line marks the position of the stress plateau which, for the present 
purpose, plays the role of a dynamic yield stress, gr^y namic ~ o.5 (this stress plateau becomes visible 

% On 



if the data are plotted in log-scale as shown in 



a linear plot, data points for shear rates 



below, say 7 = 10 4 are not distinguishable). For shear stresses below a 



dynamic 



0.5 the shear 



rates obtained from the present Poiseuille-type flow are scattered around 7 = indicating that the 



shear rate vanishes for a < a 



dynamic 



. All quantities are given in LJ units. 



Here, we study a consequence of the presence of a finite yield stress as a flow is induced in 
a planar channel via the application of an external force. The stress in such a Poiseuille-type 
flow is a linear function of the distance from the channel center. A two phase behavior may, 
therefore, occur in the glassy phase provided that the channel width is sufficiently large in 
order to ensure that stress close to the walls (where it reaches its maximum) is higher than 
the yield stress of the system. While the system response is expected to be solid-like (zero 
shear rate) in a central part of the channel (defined as a region where the stress is below the 
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FIG. 8: (Color online) Mass flow rate, Q = 2pL y Jq Z ^ 2 u(z)dz, versus applied force per particle for 
two temperatures characteristic of the normal liquid state (T = 2) and the glassy phase (T = 0.2). 
In the normal liquid state, one expects u(z) oc F e (more precisely, u(z) = F e p(L 2 /A — z 2 )/(2r])) 
and hence Q oc F e . In the glass, on the other hand, the velocity profile is blunted due to the 
presence of a yield stress. We expect the system to stop flowing for sufficiently weak external 
forces, F c < 2a y /(pL z ){~ 0.01 for the case studied here). The simulation results shown here 
conform this picture. Please note that, while the relation Q oc F e is an exact result for the normal 
liquid state, the fit Q oc (F c — 2a y /(pL z )) 2 is a purely empirical one motivated by the fact that Q 
is expected to vanish at F e = 2a y /{pL z ). All quantities are given in LJ units. 



system's yield stress), it should flow in the 'wings' delimited by this central solid-like region 
and the walls. This expectation is born out by our simulations (Fig. [3]). 

Furthermore, using the velocity gradient across the channel (Fig. E]), we define the width 
of the 'jammed' phase as the size of the region with a shear rate of j(z) = du(z)/dz < 10~ 3 . 
The accuracy of this estimate is demonstrated in the right panel of Fig. where it is shown 
that the relation W = 2cr y / (pF e ) is well satisfied by the data obtained from the above 
analysis. As to the numerical value of a y used in the above formula, our simulation results 



while the use of cr; tatlc 



0.6 



are consistent with a value of a y — (j^ynamic ^ q g (= s tress plateau upon imposed shear jsj]) 



101 ] seems to overemphasize the size of the solid-like region. 
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Our simulation results also clearly show that, as the stress increases above the dynamic 
yield stress, the contribution of a stick-slip like motion to the overall shear rate increases 
significantly (Fig. Hj). 

As a consistency check, flow curves obtained from the present Poiseuille-type simulations 
are compared to the result of previous simulations under homogeneous shear With 
the exception of relatively high driving forces where uncontrollable viscous heat significantly 
biases the present simulation results in the the vicinity of the walls, good agreement is found 
between both approaches (Figs. [6] and [7]). 

Finally, the dependence of the mass flow rate on the externally imposed force is studied. 
It is shown that the flow fully stops as F c < 2cr y /(pL z ). Again, simulated data are well 
described by this formula provided that a y = ^-dynamic « 0.5 is used. In particular, the use 
of cTy tatlc = 0.6 leads to F c = 0.0125 for a complete arrest of the flow, in contrast with our 
simulations where a finite flow is observed for this value of the external force per particle 

(Fig. ED. 

It is interesting to note that similar (blunted) shapes of the velocity profile are also 
observed in pressure driven flows of both neutrally buoyant suspensions of spheres (with 



a size of the order of 1mm) |27l. l28lj as well as red blood cells (biconcave disks of 2/im 



thickness and 8/mi diameter) |29l. |30| • In these cases, however, the profile blunting is usually 
accompanied by a migration of particles from the wall region towards the center of the 
channel ('wall migration') a phenomenon, which is absent in the case of present studies (see 
the right panel of Fig. Ej). 

Nevertheless, a qualitative similarity to the case of present simulations may be found 
when red blood cells are concerned. This similarity rests upon the fact that profile bunting 
in red blood cells occurs only if a certain amount of aggregation among red blood cells is 
present (see e.g. Fig. 7 in reference 30|). The aggregation gives rise to a finite yield stress, 
an important feature whose effects on the flow profile are the focus of the present work. 
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